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1 Introduction 



There are many phenomena that have: 1) bounds on their characteristics, 2) large number of rela- 
tively small factors composing the value of these characteristics, 3) the characteristic development 
rate is proportional to their values. To describe these phenomena, the main features of the function 
under consideration should be the following: 1) the domain of definition is the half-line, 2) the peak 
region is similar to the peak of the normal distribution and 3) the tail falls exponentially. 

These kind phenomena are in high energy cascades. To investigate the giant energy cascades 
(about ~ 4 • 10^^ eV - the conditional threshold of the relic cut-off [1]), the approximation of the 
total number of particles of the individual cascades, N{t), at depth t of the atmosphere, the function 
N^exp{—{t — / {a(t — tm) + 26^)) has been successfully used [2]. The parameters of the function 
have an obvious intuitive interpretation in this case: {tm, Nm) is the cascade maximum position, 
a is the asymmetry, h is the breadth. The physical phenomenon can be seen to have the above 
mentioned three properties as definition properties as docs the functional form: 

fen{x) ^ A{a,h)exp{-{x-cf /{a{x-c) + 2h^)) a > 0, x > c-2h'^/a or a < 0, x < c-2h'^/a. (1) 

The giant or ultra high energy cascades (UHE 4^ E'o > 10^^ eV) are rare events: 1 event per 

1 km^ sr in a year at £'o > 5 • 10^^ eV; 1 event per 1 km^ sr in a century at £'o > 1 • 10^° 
eV. Each event is examined individually. The total number of the detected cascades of energy 

> 5 ■ lO^^eV is about 15000 and those of energy £"0 > 1 ■ 10^° eV is about 200 to date (end 2007 
yr.). The must fruitful experiment with the integrated exposure mounts up to about 5165 km^ 
sr yr (Pierre Auger Observatory 01.01 2004 ^ 28.02 2007) has 11853 events Eq > ?, ■ 10^^ eV and 
only 2 events £'0 > 10^°eV. The situation due to not technical specifies experimental features only, 
but shows fiux suppression at the highest energies. The statistical and systematic uncertainness 
in the energy scale are of the order of 6% and 22%, respectively [3]. Technical difficulties cause 
registration of a part of fiux and lower quality of the processing data. 

Simulation of UHE cascades is multi-step-by-stcp process requiring large resources, so it is 
carried out by a hybrid method with plenty of internal assumptions and averaging [4]. A typical 
claimed accuracy of simulation results is ~ 5% (felt by some to be an optimistic estimation). The 
smallest change in the simulation, made by the same author, changes results by approximately 
10%. At each stage of the multistage process, a method of processing both experimental and model 
data is chosen on the basis of the quahty of this stage's input data. 

Form (1) has clear exponential asymptotic, parameterization of the fen{x) is visual and has 
correct geometrical interpretation (next section (15), (16)) besides the intuitive interpretation. 
Since the form is very simple, it may be used as basic distribution to approximate poor quality 
data. Practical usage of the refined apparatus of distribution functions to describe the phenomenon 
requires paying attention to the initial data quality. Descriptive ability of fen{x) is rather wide and 
generalizes the abilities of exponential and normal distribution (see fig. 2), therefore fen{x) may 
has wide applications. UHE cascades maximum position generated by protons in atmosphere was 
investigated using the fen{x) as probability density [5]. 

Properties of the asymmetry and of the transformation to normal distribution give ability to 
treat fen{x) for approximation of sum of finite number, n, random values with accuracy 1/n (see 
next section (22) - (26)). On analogy of the normal distribution, fen{x) can be interpreted as 
distribution of the measurement error at finite number of factors, n, which affect on asymmetric 
elementary measuring or as approximation sum of finite number of asymmetric random variables. 

The present paper deals with the methods of approximation by this form of data taking into 
account their completeness, quality and the type of the presentation independently of whether they 
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can have probabilistic interpretation or not. Approximation of any distribution by formulae simple 
and comfortable for both analytical and numerical methods is an ordinary task, often appearing 
when simulating and processing data. The exponential and normal distributions are often used not 
only because they represent distributions resulting from a physical process but also because they 
are, in our opinion, comfortable to use. One glance at a graph or histogram is enough to estimate 
their distribution parameters. Form (1) inherits the suitability. 

Experimental or simulation results are usually represented by the mean value and variance only. 
Often, however, the data, considered as random variables are asymmetric. Therefore, representation 
of data by the mean value and variance only is insufficient. The next step is representation of the 
variance as a sum of left and right parts, and form (1) is convenient for the representation. The 
present article intends to show that the asymmetry can be taken into account even in a case of 
poor quality data. 

At the beginning of the article, the properties of form (1) are presented as the simplest case 
between the exponential and normal distributions. Then methods of approximation are developed 
with taking into account the initial data quality and, finally, the methods are used to study UHE 
cascades. These last sections (thought may be interesting in itself) illustrate use of the methods. 
In our opinion, the elaboration has general character and may be useful in other fields. 



2 Some properties of form A{a,b)exp{—{x — cY/{a{x — c)-\-2b'^)) 

The parameterization is convenient for interpretation of the parameters. For standard support 
X G (0, cxd) the form is B{s,q)exp{—{sx + q/x)/2), where s > 0, q > 0. The distribution fen{x) 
is infinitely-divisible [6, b] and belongs to the family of Generalized inverse Gaussian distribution 
(investigated by B. Jorgensen [7]). The probabihty density of the family is 

where a; > 0, ATp is a modified Bessel function of the third kind, with index p < ±oo, s > 0, and 
q>0. 

The distribution fen{x) is a special case [6, a] of Hyperbolic distributions: its hyperbola has 

vertical asymptote. The special case generates the family Hyperbolic distributions (introduced in 
and investigated by O.E. Barndorff- Nielsen [6]) by normal variance-mixing. The probability density 
of the family is: 

f (x) = ^-^/({x-^)^+q){0^+s)+{x-^)0 

K,{qs) 

where x e (— oo,+oo), s > and q > are parameters. The family is normal variance-mean 
mixture Y — + PV + ay/VX at cr = 1, and P are free parameters. Here random variables 
X and V are independent, X is normal distributed with mean zero and variance one, and V has 
probability density fen{x) on support x G (0, oo). 

Formulae for the general characteristics (normalization, moments, Laplace transform, Fourier 
transform) of the distribution fen{x) may be obtained with help of formulae 3.471.9 of Gradsteyn 
and Ryzhik [8]: 

POO 

/ xP-^exp{-{sx + q/x)/2)dx = 2{s/qY/^Kp{^) Re{s) > Re{q) > 0, (2) 
Jo 

where Kp is a modified Bessel function of the third kind, with index p. 
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In particular, for normalization factor of fen{x) we have 

A{a,b) = [a{2b/afexp{{2b/af)Ki{{2b/af)]-\ (3) 
Asymptotic expansions [8, 9] for Ki{z) 

Ki{z) ^ - + ^In^ + o{z^) z^O, 
z z z 

if.W = y^e-'(l + l-j^+...) .^oo (4) 

show that the feni^) tends to the exponential distribution, fe{x), as 2; — >■ and to the normal one, 
fni^), as 2; — > 00. distributions. 

At scaling x ^ x/a with saving equality -f{-,a) = f{x,a') the parameters of the fe{x) and 
fn{x) are transformed a' = a ■ a. The parameters a and b of the fen{x) are transformed in the 
same way. Thus, a/b = a' /b' therefore this ratio is scaling invariant and can be considered as the 
essential parameter of fen{x)- It is convenient to define the scaling- invariant parameter 

-6*V. (5) 



take it as a parameter of type of the distribution feni^) a-nd write the normalizing factor (3) in the 
form: 

A{a, b) — [a ■ z ■ exp{z)Ki{z)]~^ . 

Obviously, the parameter z may be treated as a measure of asymmetry (skewness) of the fen{x), 
but at data processing it is convenient to use another scaling-invariant function also. 

We define af and cr^ as left and right parts of the variance cr^ (above the mean value) by the 
formulae: 

/•oo 

xf{x)dx, (6) 

ic-262/a 

+00 X +00 

a^^ j{X- xff{x)dx = af + ^ j {X - xff{x)dx + j - xff{x)dx. (7) 

—00 —00 X 

The ratio 

= a^/al (8) 

is invariant under scaling of x and therefore depends on z only, and also may be taken for 
measuring of asymmetry. For exponential distribution R^ = e/2 — 1. Since R^ for /en (2^) is 
monotonic function bounded between its values for fe{x) and fn{x), we have 

e/2 - 1 = 0.359.. < i?^ < 1, (9) 

for a > 0. Fig. 1 shows the measure of asymmetry, 1 — crf/a^, of the fen{x) from exponent to 
normal distributions. 

One can sec that fen{x) is a unimodal, bell-shaped, right-asymmetric (a > 0) distribution. To 
illustrate the descriptive ability of fen{x), fig. 2 shows the distributions at c=0, 6=1 and variable 
a=0, 1, 5, 10. One can see that fen{x) changes from fn{x) towards fe{x). The dotted line in fig. 2 
is the easily recognized standard form of the normal distribution at a=0. An examination of the 
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figure shows the distributions at a 7^ have supports on half-hnes, with the left ends located on the 
hyperbola a = —2/x according to the conditions in (1). The maximum density (equal to the value 
of normalizing factor A{a, b)) decreases with increasing a. For a =5 the value A(5,l) is singled out 
in the figure. The distance, x, between the mode, c, and the mean value, X, 

X^x + c (10) 

increases with a gradually because the distributions are steepening on the left side and flattening 
in the right side. It can easily be imagined how the form of the distribution is changed during 
scaling and see that no one of its can be obtained from other by scaling. 

Thus, probability density fen{x) defines a one-parameter set of distributions with rather wide 
descriptive ability. 

Note that the mode is a more natural characteristic of the distribution feni^) than the mathe- 
matical expectation, whose distance from the mode characterizes asymmetry. Such kind of asym- 
metry is known as Pearson's skewness: {X — c)/a, which may be obtained from lines 3 and 4 of 
table 1 (see below). 

It is interesting to study the width of (1), W, with its left part, Wi, and right part, Wr, at the 
mode: W — Wi + Wr (see fig. 3). The width of distribution is used quite often. For example, 
an error bar graph is usually used to characterize the accuracy of measurement results. The —a 
and +a values represent a normal distribution width at the e~^/^ level, i.e. at the level — e^/^ 
where the probability density is less by a factor of e^/^ than at its maximum. In some cases the 
width of distribution has a physical interpretation and measures a physical quantity, for example, 
a Cherenkov light lateral distribution, or the lifetime of a resonance. The width is less fluctuating 
characteristic then other and is used as robust characteristic. For UHE events heavy nuclei are 
expected to produce smaller shower- to- shower fluctuations then protons. Therefore the width of 
distributions of its characteristic is sensitive parameter. These properties are useful and are used 
wide [10.11.12]. 

These W, Wi, Wr at the level d is defined by solutions of the equation: 

A{a, h)exp{-x^ / {ax + 26^)) = A{a, h)/d = A{a, b)exp{- In d). 

One can see from the equation x'^ — x • alnd — 2b'^lad — (see also fig. 3) that 

1. The width of the distribution is determined by the discriminant of the equation: 

W = ^a^{lndy + 8b^ Ind. (11) 

2. The parameter a is determined by the difference of the right and the left parts of the 
distribution width (at the mode c — 0) i.e. 

(Wr-Wi)/\nd. (12) 

3. The parameter b is determined by the product of the left and right parts of the width 

262 = WiWr/lnd or b = (1^,1^^/(2 In d))^/^. (13) 
For future convenience rewrite (11)-(13) for the level d = e: 

W = Va2 + 862, (14) 

a = Wr - Wi, (15) 
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= WiWr or 6 = {WiWr/2y/^. (16) 

Formulae (15), (16) have evident geometrical interpretations (see also fig. 3) which induce a 
physical interpretation of parameters in particular cases. 

Note the relationship between the widths and the parameter z (5) defining the distribution form 



Obviously, that the width at one level is not sufficient to determine fen{x), but adding the ratio 
parts of the width (divided by the mode or the mean) is sufficient, as it is shown by (11) - (13). 

One can see some analogy of ways to describe the distribution function in terms of width and 
variance, because both of it's describe a breadth, but width is function of the level d. It is possible 
to create approximation of UHE cascade profile in term of width [13]. Mistakes caused by the habit 
of believing that the width of distribution is defined by the standard deviation should be avoided. 
For example, if it is assumed that W — 2a (at the level d~^ — e"^/^ ), i.e. if the relationship 
between the width and the variance is assumed to be the same as for the normal distribution, one 
obtains a formula for the variance: cr^ = 6^ + 0^/16, that is true at a = only. 

The same as for normal distribution often it is needed cumulative form of the fen{x) distribution. 

The function 

/-^ exp{-{x - cf / {a{x - c) + 2}?)) ^ 
^en[X,a,b,c) - J_^^^^^^^ a{2h/ayexp{{2b/af)m2h/ayr ^''^ 

is the integral form of one-parametrical distributions fen{^)- Since, for practical using, it is conve- 
nient to have (17) together with the table binding the values of parameters a, 6, c with the basic 
characteristics of the distribution (the normalizing factor, the mathematical expectation, the vari- 
ance, and skewness ), the standard form fen{^) is taken in form: 

F.„iy,.)^r ff .. (18) 

J-z/2 z ■ exp{z) ■ Ki{z) 

corresponding to a = 1, c = 0. Thus, 

a, b, c) = Fen{{X T c)/(±a), {2b/af) = Fe„(y, z), (19) 

where the sign at the parameter c is opposite to the sign of the parameter a value. 

So, for one-parametrical distribution feni^) we have table 1 with lines: z, A\a\, x/a, a^jc? — 
1, (yf/crl, y, where the vector designates the lines Fen{y,z) (18) for the y values, listed in the 
first column. After calculating z = {2b/ a^; the values of the normalizing factor, mathematical 
expectation, variance, and the ratio of components of the variance can be obtained using the values 
of parameters (a, b, c) from the top 5 lines; and conversely, the values of mean, the variance and the 
ratio of components of variance can be used to obtain the values of ^4(0, b),a,b,c (see below in section 
3 the description of method A and the numerical example in section 5, footnote 3). Dealing with the 
cumulative form (the distribution function) is defined by formula (19) determining Fe„(X, a, 6, c) : 
it is necessary to find in the table the value of Fen{y,z), where y = {X^c)/{±a),z = {2b/ a)^. The 
numerical example see in section 5, footnote 4. 

Let us present the list of the basic properties of the fen{x) as bridge between the exponential 
and the normal distributions. 
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Distributions 

a w 



exp{-{x-c))/a) exp{~{x-c)'^ /{a{x-c)+2b'^)) exp{-{x-c)^ / {2a'^)) 

{2b/afexp({2b/ay^)Ki{(2b/a)'^) 

support c < X < oo —2}? I a + c<x<oo oo<a:< +00 

Symmetry - - + 

Interpret, of parameters +, + +, +, + +) + 

Scaling + +, + + 

Unsealing parameter - hja = (26/a)^] - 

Mode c c c 

Characteristic Function ^ ^'''^ e-^+^*^ 
Mathematical 

expectation X c + a c+ —{K2{z)/ Ki{z) — 1) c 

Variance = {X - cf 462[(||i)2(_|) + + |] ^2 ^^2 

moment /.3 20^ a462[(gM)3(|) + (gM)2(_ 3^^) + M£l(_^ + 3) + f^^] 

Ratio crf/al, (a > 0) e/2 - 1 e/2 - 1 < (Ji/cr^ < 1 1 

Width at level e-^ a ^/a? + 86^ [W/,^ = if a/2] 2^2(7 

One can see how properties of the normal or the exponential distributions for the intermediate 
case under consideration are generalized, transforming to appropriate relations at a — > or 6 — > 0. 
The fen{x) enables description of some asymmetric distributions. The possibility is furnished by 
parameter a but, in general, it is incorrect to think of this as a skewness parameter, because its 
value is not scaling- invariant. Any monotonic function depending on a/6 only may be taken as a 
measure of asymmetry (skewness) for fen{x), for example 1/z defined by (5) or 1 — i?^ defined by 
(8). 

Note one more property of fen{x) caused by transformation of feni^) to normal distribution. 

Since under weak conditions, sum of n symmetric random variables is approximated by the 
normal distribution with an accuracy ~ l/n, but asymmetric ones with accuracy ~ 1/a/^, that it 
is possible to find such asymmetric distribution which approximates the asymmetric sum with an 
accuracy ~ 1/n. One of the such kind distribution is fen{x). 

Indeed, let /(x)s be the probability density of x = Y17=i where are the independent 
random variables with a common distribution which mean value is m and standard is a. Write 
Edgeworth expansion [14] for the f{x)-E including term R^n — 0{n~^): 

f{xh = ^{x) - n-'^^^^,^^'\x) + (^72^('H^) + ^7iV^')(^))(n-^) - (20) 

2 

where (p^''\x) is k-th derivative of >^{x) = ^=e~^, 71 = /^3/cr^ and 72 = /i4/o"^ — 3 are the skewness 



/2n 

and excess of the distribution. Reminder that ip^''\x) = {—!)'' Hk{x)ip{x), and Hk{x) is k degree 
Hermit polynomial. Using fen{x) for approximation of f{x)-£, we have the discrepancy 

A = ip{x) - n- 5 17,(^(3) (^) _ A{a, b)exp{-{x - cf/{a{x - c) + 26^)) + 0{n-^). (21) 

We will take such parameters a, h, c that term with will be covered by fen{x). From (3) and 
(4) we have A{a, b) = -^^^ + O(a^). Dividing (21) by (p{x) and using series for obtained exponent 
one can has: 

(A - 0{n'))l^{x) = 1 + n-^i7i(^^ - 3x) - ^(1 - ^^j^^^f^^^ + y ) + 0{a^) + 0{c') 
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il , o , 1 2x^ - 4:xc + 2c^ - x'^a(x - c) - x^2b'^ ^, ^,9, 
= "''S^-'^ - - 6 + 26(a(.-c) + 26^) + ' + > 

_i 1 ,0 „ , ax^ Acx 

— n 2 — 7i(a; — Jx) — 



3! " ' ' 2h{a{x - c) + 262) 26(a(a; - c) + 2b^) 

(262 - 2 - ac)x^ , 2c2 1 , . , 2^ , 2, 



2Ka{x - c) + 262) + 26(a(.-c) + 262) " 6 + ^ + ^^"^ + ^^^^^ ^'^^ 
The equation shows the needed values of the parameters a, 6, c: 

2 1 ill 

a= -7in-5 + 0(n-^), 6 = (1 + 0(n"i)i, c= — 7in"^ + 0(n-^). (23) 

3 2 

Then by direct substitution one can be convinced that A/(p{x) — 0{n~^). That is 

f{xh = fen{x) +A + R4n = fen{x) + ^{x)0{n-^) + i?4n = fen{x) + 0(71"^). 

According to (23), the parameters a, 6, c become to depend on n^^^'^ji, therefore the support 
(1) of the fen{x) depends on n~^/^7i: x > c — 2b'^ /a = — 3n^/2^-i _ rr'^l'^^i\l2. Out of the region 
we have: < /(a:)s < 0(n~^) and we define fen{x) = 0. The problem because of negative values 
of H^ix) [14, 1946; 6, 1989] does not arise. 

Thus, the fen{x) uniformly approximates f{x)s for — oo < x < oo with accuracy 0(n~^). 

If have not common distribution [14, 1937], the 71 and 72 must be chosen by 

i—\ i— 1 i— 1 i— 1 i—1 

where and are the third and the forth central moments of ^j, respectively. 

The property has various consequences and seems to be useful for practical use. 

As one can see, (23) can be defined more exactly. For example, one can obtain the fen{x), which 
three first moments coincide with the same of the f{x)s- 0, 1, n'^^^^i. Prom (2) and (4) one 
can obtain the moments of the fen{x) in form: 

, , , 3 3 (3/ 3 (3/ 

Mn{x)) = ^06^ + ^a=^ + ill . ^ + . . . , 

and then the parameters of the fen{x) in form: 

2 1 2 3 5 
a = gTi^"^ - ^7?^"^ + C>(7i^n-i), 

&^ = l-^7^- + ^7^- + 0(7M, 

1 _i 3 Q _3 c _5, 

c^-^Ti^ '^2332^^1^ 2+0(7in 2), 

when n^^'ji < 1. 
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For better accuracy of approximation of f{x)j: by the fen{x), in specific cases we can introduce 
free parameters ai,bi,ci: 

a = -9 + ai9^, b'^ = l + bi9^, c=--9 + ci9^, 
3 2 

where 9 = n"^/^7i < 1. and take into account the last term of Edgeworth expansion (20), which 
is ~ 0{n^^). As one can see, it is impossible to eliminate R^ni^) in general case, but the line 
{\ai — |6i — Cijx + which is part of the i?4„(x), can decrease -R4n(x) ~ 0{n^^) in any region. 

The property (approximation of finite number of asymmetric random variables) allows to inter- 
pret the fen{x) analogously with the variant of the popular interpretation of the normal distribution 
as measurement error distribution. ^ is elementary measuring (^ — m is elementary error). Asym- 
metry of the ^ may be caused by either object of the measuring or by device of the measuring. If in 
practice the finite number, n, of essential factors (but not infinity number of infinitesimal) affects 
on result of the measuring, x, that asymmetry of the ^ is not terminated and the distribution 
of X is asymmetric. The fen{x) approximate this asymmetric distribution with accuracy ~ 1/n 
because can not take into account many other characteristics of the ^. (Infinite n terminates these 
characteristics together with the skewness)^. 

It is obvious, that not only the fen{x) has this property. Interpretation of the parameters and 
simplicity of use of the normal and exponential distributions is inherited by fen{x). 

For illustration, we shall examine a concrete example of the study of the cascade maximum 
depths in the atmosphere using the information from [15]. P{tm) is the distribution of the depth 
of maximum, tm, of 500 simulated cascades generated by protons with energy lO^'^ eV at depth t 
(g/cm^) in the atmosphere, shown as the histogram presented in fig. 3. 

One can draw by hand a smooth curve approximating the histogram (see the dashed curve in 
fig. 3). It is a right-asymmetric unimodal distribution. Let us drop perpendicular from the curve 
maximum (670, 35) to the abscissa and obtain the value of the mode c = QlQg/crn?. Then draw 
the horizontal line approximately at the level 1/3 (e~3) of the maximum of the curve and mark 
the values of Wi=bS and Wr=M. Then (15, 16) yield a=94-58=36, h = (58-94/2)1/2 ^ 53 2; so the 
analytical approximation by form (1) can be written: 35exp{—{t — 670)2/(36(t — 670) + 2 ■ 52. 2^)). 

The distribution does not contradict to the physical understanding of the problem and the 
complete description in [15]. One can test accuracy and acceptability of the formula by the same 
way as for using the normal distribution. For example, it is very simple to compare the mean value 
of the approximation with that calculated directly, 699.63 g/cm^. For this, one can calculate the z 
value defined by (5): z = {2- 52.2/36)^ = 8.41. Then, using the values of z and a, from table 1, one 
can obtain x = 0.733 ■ 36 = 27.83; so account of shift (c=670) in (10) yields the mean value of the 
approximation: X = 670 + 27.83 = 697. Sg/cm"^. (Such good fit for a hand drawing is accidental). 
To obtain the probability density function P{tm), it is only necessary to know the normalization 
factor. It can be obtained from table 1 by using z and a: A{a, h) = 0.260/36 = 7.2 ■ 10^^. 

This example illustrates a qualitative analysis and using the table 1. The distribution may be 
used as the first approximation for the correcting process. Regular methods considering specifics 
of the initial data are described below. The efficiency of application of particular properties of 
the distribution depends on many things: on the available set of the data under conditions of a 
concrete problem and even on the form of representation of this information. 

^Authors can't point out physical process which the fen{x) describes exactly. 
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3 Methods of creation of approximation 



Although it is obvious that fen{x) can be determined by the first three moments and the general 
methods are known, we focus our attention to the problems with the incomplete or poor-quality 
data. It is a typical case in the investigation of complex objects, for example for UHE cascades. 
Since the cascades are rare, difficulty registered and one can see only part of the cascade (in 
particular, it is longer than deep of atmosphere), input data may have errors about 30% and may 
be incomplete. For example, some details may be seen in [11, 16]. Exclusions are so-called 'golden 
events' which are used as etalons. Thus, for accounting asymmetry use a third moment may be 
incorrect (or in any case is unreliable). The above properties of fen{x) make it feasible to create 
regular methods of approximation (in contrast to the by-hand method described in the previous 
section) accounting for asymmetry without using third or higher moments. 

Method A. Approximation with conservation of mean, variance and ratio of the variance's parts. 

If the mean and variance are calculated, obtaining the parameters of the exponential and normal 
distributions is reduced to the comparison of the calculated variance with the variance of a standard 
form (for determination of the scaling parameter) and to the comparison of the calculated mean 
with the mean value of the standard form (for determination of the parameter shift). Knowledge 
of the values of the mean and variance is not enough to determine the parameters of fen{x), since 
the standard form of this distribution is not unique but is one-par ametrical series depending on a 
scaling-invariant parameter, which defines type of the distribution and may be taken as measure of 
asymmetry. Therefore, if the value of a scaling-invariant parameter is found first, the problem of 
obtaining the values of the parameters of fen{x) will be reduced, in essence, to the same well-known 
procedures that are used when dealing with the well-known neighboring distributions. 

In this connection, one can use the measure of asymmetry introduced by (8), since the value of 

can be found using the initial data. This value allows one to find (for example, with the help of 
fig. 1 or table 1) the parameter z of the distribution, i.e. to solve the problem, remaining within 
the framework of habitual procedures of work with the well-known distributions. 

Formula (8) specifies an intuitive estimate of (the so-called method of substitution). Tak- 
ing into account continuity of all functions under consideration, the estimate will obviously be 
consistent. Its properties will need to be studied in more detail before wide practical use of fen{x). 

So, using the initial data, 3 numbers are calculated: the mean value, X, the left and the right 
parts of the variance, af, a^. Then value of R^ (8) determines the value of parameter z with the 
help of table 1. The value of z determines both values of variance, [cr^/a^], and mathematical 
expectation, [x/a], of the standard form of fen{x), [a = 1, c = 0]. A comparison of the calculated 
variance with its obtained value, [a^/a^], determines the parameter a = a/ {cf + a^) /[a'^ / a?]. Then 
similarly by taking into account shift, we obtain: c = x — \xlo\^ {af + cTj.)/ [ct^/ a^] . To complete 
the discussion, the parameter b is obtained from (5); and the parameter A is obtained from (3) or 
with the help of table 1. ^ |^ The distribution fen{x) with the parameters calculated in this way will 
be constructed to have the mean, variance and ratio of the variance components of the data. (The 
numerical example is presented in footnote 3 of section 5). 

This method is suitable for various presentations of the initial data; for example, as histograms 
(interpreted as density of the distribution function, or without such interpretation and considered 
as values of some function obtained by simulation of a phenomenon), or as a set of results of direct 
measurements of some physical characteristics; since the method only requires it be possible to 

^ These 3 numbers (X, af, cr^) are representation of fenix). The numbers (X, ai,ar) are convenient for graphical 
representation of the asymmetrical distribution function, which is similar to traditional representation results of a 
measuring by mean value and error bars. 



10 



calculate the mean and the left and right components of the variance. For direct measurements the 
stage of creating histograms can be omitted in data processing. The method can be used for poor 
statistics (say less then 50) when it is difficult to obtain a reliable histogram. The procedure of the 
method can be easily inserted in a popular procedure of previous data processing (calculation of 
the mean and dispersion). At last, the method allows checking the applicability of (1) by testing 
the equation (9) for the asymmetry. 

Method B. Approximation by two pairs of parameters. 

Let us study a way of approximation based on a specific property of form (1) using the linear 
(at x) form of the denominator in the exponent: 

a{x - c) + 2b^ = (x - c)yin{A/f{x)). (24) 

It can be seen that the approximation f{x) (in particular, the histogram), at the known position 
of the maximum of (1), i.e. at the known value of pair parameters A and c, is reduced to the linear 
approximation of the computable right part of (24). 

The pair {A, c) determines the position of a curve (1) on the coordinate plane and its central part, 
but actually the form and scaling is determined by the pair (a, 6) (see (5), (15), (16)). The initial 
data and their processing often differ in the central and in the peripheral parts of the distribution, 
because different processes dominate in the central part and at the periphery. Both for physical 
and statistical reasons, the central part is usually less subject to fluctuations than the periphery. 
In the case of small fluctuations in the central part, it is possible to obtain enough accuracy the 
values of A and c by a fltting maximum of the region (e~^ = 1 — x'^ + x^/2 — ...) to a parabola 
which approximate, say, 5 points at the region of maximum. It is possible (and frequently useful) 
to exclude some of these points from the subsequent calculation (24) when obtaining a and b. An 
approximation constructed in this way is based on a concept of a mode (instead of the mean as in 
method A) as a more natural and robust characteristic. 

Fig. 4 illustrates method B for a simulated individual shower. The shower was generated by a 
proton with energy lO^^eV and zenith angle 9 — 44.4°. The circles show the input values N{t), the 
sohd line is for approximation (1) and the squares are for the values of the right part of (24). The 
area near the maximum is presented for more detail in the inset. The parabola determining tm and 
N(tm) is given by the dashed line. For the linear (at t) approximation determining the values of a 
and b 

a{t - tra) + 26' = {t- t„,f/ln{N{tm)/N{t)) 

the ordinate axis with natural more detailed (~ 4 order) scale is located to the right. From the 
figure it can be obtained that Nm = 7.39 • 10^ tm = 740; a a; (1.43 • 10^ - 7.8 • 10^)/1000 = 
65, 26' ^ 1.43 • 10^ - (1400 - 740) = 1 • 10^ So we have: 

N{t) = 7.39 • W^exp{-{t - 740)V(65(t - 740) + 1 • 10^)). 

This method is suitable for truncated data, i.e. in the cases when the tail (or nose) of a distribution 
is unknown, or for incomplete data with "holes" (for example, when an intermediate part of whole 
diapason is inaccessible for measuring or a histogram columns are lost during data transmission) 
since the method requires calculation of the maximum and the linear approximation. From fig. 4 
it can be easily estimated how the removal of some points will affect the approximation, so one can 
judge the quality of the approximation. 

Method B is fitting method essentially and based on the robust characteristics of the data. For 
the case of sufficient statistics the method has an advantage over method A, which can be used for 
complete data only. 
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Concrete conditions may lead to construction of methods intermediate between methods A and 
B. For example, if the mode and mean can be determined quite reliably while the parts of variance 
seems unreliable, the K. Pearson skewness {x — c)/(t [14, 1946] based on the mode can be used to 
construct fen{x). 

Results of the methods can be used as starting point for refinement by fitting procedures. 

4 Some general problems application for UHE cascades 

Ultra-high-energy (UHE) events are very rare, so this fact leads to two consequences: a limitation 
of using the usual probability theory apparatus and the highly important role of the Monte-Carlo 
method. Since the phenomenon is complex, the central processor unit time required would rise 
in proportional to primary energy and a single shower could take well over a century at 10^° 
eV ([4] 1997). Practical modehng is carried out by a hybrid method that introduces unphysical 
fluctuations. Therefore when analyzing the data, it is important to pay attention to the quality 
input data for each step of the multi-step data processing. In this section we discuss in brief some 
general problems of the application of form (1) for the UHE cascades. (For some more details see 



The results of Monte Carlo simulation of air showers generated by protons with energies of 
10^^, 10^° and 10^^ eV for zenith angles ^=0°, 24.6°, 44.6°, 60° are used. Calculations were carried 
out in the framework of the quark-gluon-string model [17]. The hybrid method [18] was used to 
estimate both the mean values and the standard deviations of the parameters under consideration. 
The method enables accounting for fiuctuations of the number and the location of interactions in 
the atmosphere and energy release of primary protons, using the Monte-Carlo procedure, while the 
development of cascades from numerous charged pions is considered in an average of step-by-step 
approach [19]. For energies higher than the threshold value of 10^^ - 10^'' eV (depending on the 
position in the atmosphere), the propagation of electrons and photons was simulated with account 
for Landau-Pomeranchuk-Migdal (LPM) effect suppressing of bremsstrahlung and pair production 
[20]. (Details of the simulation for Yakutsk experiment can be seen in [21]). The results obtained 
were presented as functions at the atmospheric depth t with step At =50 g/cm^. The technical 
quality of the simulation results was rather high; local unevenness and unphysical fluctuations are 
rare, their deviations rule, less than the natural cascade fluctuations. 

For approximation of UHE cascades use is made of various formulae. A well-known Gaisser- 
Hillas formula was introduced by [22] to describe the mean particle number N{EQ,t) in a vertical 
{9 — 0°) shower generated by a proton of energy Eq > 10^^ eV at depth to: 



(t' is measured in units A=70 g/cm^ from the generation depth to; the depth of the maximum 
is t^=0.51 ln{Eo/c) — 1 and Nq=0.0A5, c=0.074 GeV). This formula was used later to approxi- 
mate individual cascades, assuming it to comprise four parameters: tm, for the position of the 
maximum and to, A are free parameters. 

Other functional forms have also been used; e.g., two parts of Gaussian functions with two 
different widths ai and a2 is used for t < tm and t > tm, respectively, [23]; or other forms [16], for 
example, above mentioned parameterization using the width of the distribution [13]. In all cases 
the approximation accuracy should be much less than physical fluctuations, since frequently the 
fluctuations are the object of interest when UHE events are investigated. 



[5]). 




(25) 
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We use the approximation of giant cascades developed in [2]. The number of each sort of 
particles (electrons, muons and their sums) in each individual cascade was approximated by form 
(1) as a function of t > 0.5t„t with an accuracy about 1.5%. In particular, for the total number, 

(t) , of particles at depth t in an individual cascade number i 

N'{t) = N^^exp{-{t - ej/{a\t - fj + {2h'f)) at t > 0.5C, ^ = 1 - 500, (26) 

where t^, N^f^) are the maximum of the cascade, a* is the parameter of asymmetry, and ¥ can be 
thought a breadth parameter. Approximation (26) has been shown in [2] to be stable. The form 
describes well both the mean and strongly fluctuating individual cascades. 

We describe each individual cascade by 4 parameters; an individual cascade is a concrete re- 
alization of 4 random variables with probability density denoted as P(tm), P{Nm), P{a) and P{b) 
which describe the entire ensemble of cascades. This approach leads to the definition of a mean 
cascade as a cascade with parameters each of them mathematical expectation. Such definition is a 
natural generalization of the traditional definition of a mean cascade where the number of particles 
at the depth t takes one its expectation value. 

The function P{t„i) shows a narrow peak at the left side from the mean value, i„i, and a 
wide range (~ 6a) of the parameter with approximately exponential asymptotic attenuation. The 
function P{N^ has the parameter range ~ 5(7, the mean shifted to the right and the maximum 
to the right from the mean, with a peak wider than that of P{tm)- The function P{a) has the 
parameter range ~ 4(j, its form is similar to P{Nm)- The function P{b) is similar to P{tm), the 
mean value and the maximum are shifted to the left. The distribution functions describe the model 
results in full detail in a compact form (for example see fig. 6 below). The functions are very 
sensitive: small changes of the simulation (whose influence is difficult to see) are noticeable in the 
distribution functions. 

Since cascade processes dominate in the UHE cascades, asymmetry of the cascades is not so 
large {z >> 4, see for the example fig. 4), that N^{t)dt ~ \/2TTN^b^ , see [2] for details. It is 
useful to define a new parameter s* = \/27riV^6^, which fiuctuates less than and may be used for 
primary energy determination. Analysis of the distribution functions reveals showers with super- 
fluous characteristics. Special investigation of RiO.5% outlying showers shows some shortcomings 
of the model and simulation techniques and a way to correct them. The parameter s is the most 
useful for revealing the unphysical fluctuations. Thus, the known interpretation of the parameters 
enables checking the simulation. 

In cases high quality of the data makes it feasible to apply both A and B methods. However, for 
higher energies and smaller angles, as well as for consideration of some characteristics of cascades 
(for example, muon flux) the data will be incomplete. Vertical cascades of higher energies poorly 
reach a maximum at the sea level, so only half of the distribution can be observed. Besides that 
it is necessary to take into account restricting outer edges of (26) (t > 0.5tm) and the increase 
fluctuations at the cascade's peripheries. When both methods are applicable, they give practically 
the same results: the difference in tm is, as a rule, less than 3g/cm^. (For comparison, accuracy of 
the parameter at experimental measuring is 10 g/cm^ in the best special cases [11]). 

5 Distribution function of the depth of cascade maximum 

The shower maximum position, tjn,Nm, is measured directly [11]. The characteristic is used for 
estimation of primary energy and mass composition of primary cosmic ray flux. The cosmic ray 
composition is studied by comparing the observed values with predications from simulations for 
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different nuclei. The change of the observed tm distribution is used to derive estimates of the change 
in primary composition. At a certain energy the average value tm of primary mass A, is related to 
mean im of the proton. 

Physical basis for such way investigation is domination of proton component (~ 80%) in flux of 
cosmic rays and known differences between cascades generated by nuclei and photons and cascades 
generated by protons. Then briefly speaking, due to larger multiplicity in firsts interactions of 
the showers generated by nuclei develop faster and less fluctuate than proton ones. Therefore 
(at the same energy) the distribution function of tm from nuclear primary is shifted at less depth 
(~ 7Qg/crn?) and has less width (as think the authors, the skewness must be less too). But cascades 
from photon primary due to less multiplicity difference from proton ones to the opposite side, i.e. 
its depths are larger and distribution functions are wider. 

Thus, tm of proton generated showers is basically characteristic for energy and mass composition 
estimates. For these estimations are used the showers which estimations of uncertainness tm and 
total energy are smaller then AQg/cm? and 20%, respectively [12]. 

The depth of cascade maximum is determined by two physical processes: by the depth of the 
first interaction (with a known exponential distribution) and by fluctuating cascade development 
(from the first interaction point) which resembles the central part of the normal distribution. The 
distribution function of the depth is not a convolution of the exponential and normal distributions, 
but rather the function must be determined on the half-line, has a maximum region similar to the 
normal distribution, and has an exponential attenuation. Therefore, it is natural to use fen{x) as 
an approximation of the probability density of the cascade maximum depth, P{tm)- As will be seen 
below, the first process dominate since the asymmetry is significant and invariant parameter of the 
approximation z < A. 

To estimate P{tm), we use described in previous section simulated cascades of energies of 10^^, 
10^° and 10^^ eV for zenith angles ^=0°, 24.6°, 44.4°, 60°. The cascades are first approximated 
by fen{x) using method B. This representation yields tm for each of the simulated cascades and 
the total set of {tm} is input data for the step of obtaining a distribution function. Note for the 
second step the quality of the data may be worse than input data for the previous step, because an 
error of each tm can lead to unphysical fluctuations: local groups and rarefactions. The statistics 
may be poor (especially for natural events), but the data have not any 'holes' and cut of the tail. 
Therefore for the step we choose method A just to the initial set {tm}- (As recommended by section 
4, histogram of {tm} is not required) 

For the case Ep = lO^^eV, ^=44.4°, the probability density P{tm) with the parameters A{a, b) = 
0.0147, a = 38.7, b = 22.1, c = 752.5 is presented in fig. 5 by solid curve. Though the distribution 
is obtained without any histogram, for visual illustration of approximation method A quality there 
is one histogram of tm which is given by solid straight-lines. Dotted histogram is experimental 
data [10] shifted on lOOg/cm'^ to the great depth. The discussion of the experimental data will 
be below. The isolated fragment is representation of the distribution P{tm) through mean value 
and components of variance. The representation was discussed in connection with method A and 
mentioned in footnote 2. Note that accuracy of mean value may be presented by diameter of the 
point which denotes the mean value. (The representation gives ability clear to see the curve of the 

^For example, take case, Ep — lO^^eV with zenith angle 6 = 44.4°. In accordance with method A (section 3), 
calculation gives the mean value tm — 785, the left and right components of the variance: sf — 514 and = 1255. 
We calculate = sf /s'^ = 0.41 and, using fig.l or 5— >1 lines of table 1, find the value z = 1.3. For this z, table 1 
gives: A{a,b)\a\ = 0.568, x/a = 0.841, a'^/a'^ = 1.18. Then it is consistently calculated that = (514+1255)/1.18 
= 1499, (26)2 ^ 1.3/1499^ (a = 38.7, b = 22.1), A{a, b) = 0.568/38.7 = 0.0147, x = 0.841 • 38.7 = 32.5, c = 785-32.5 
= 752.5. 
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fen{x) distribution as —a ■ +a to see normal distribution). 

The construction of the probabihty density P{tm) was carried out, in essence, in an integral ap- 
proach, avoiding a histogram, therefore, it is natural to use the Kolmogorov's test. The goodness of 
fit between the empirical distribution function F*{tm) (sum polygon) and the obtained distribution 
function F{tm) 



yields: A = -/n max|F(tm) — F*{tm)\ = 0.59 and, accordingly, the probability of the observed 
discrepancy, P{\) ~0.88, may easily be explained by statistical fluctuations. A similar consideration 
was made for other energies: 1020(e1/) ^ P(A) ^ 0.71, W^{eV) =^ P{\) ^ 0.72. 

As soon as statistics (500) permits to check the A method, let us verify the results by x^- 
test. We use the histogram with columns presenting probability equal to 1/k, where k is the total 
number of the columns and let A; = 20 (example for Ep = 10^^ eV see in fig. 5). Pearson's 
test of the simple hypothesis (a = 38.7, b = 22.1, c = 752.4) against an alternative yields 
Pwix^) ~ 45% (for 19 degrees of freedom), then for composite hypothesis x^-minimum method 
yields a — 41.5, b — 17.9, c = 749.5 and Pieix^) ~ 70%. The results for simple hypothesis are 
listed in the left, but ones for the composite hypothesis are listed in the right side of table 2. When 
analyzing the table, we take into account that the results for 10^° are less accurate then others. 

Due to UHE events are rare; comparison with experimental data is difficult. We can take 
distribution of tm for natural all-energy, all-zenith angle and all-kind primaries [10], only. As it 
will be shown below (see fig. 6) we can neglect the angle dependence. Approximately 20% of the 
mixture is iron which has tm less then ones of proton 85 gcm~^. The energy dependence may 
be taken into account partly and only roughly. Energy spectrum of flux in the region of so-called 
'ankle' (~ 2 ■ lO^^eV, see fig. 7) is complex (but about E~^); in any case the majority of the 
mixture consists the cascades of low energies (~ 10^^ eV). Increasing average tm with increasing 
energy at one order (so-called 'elongation rate') is about 71 gcm"'^ / decade at 10^^ eV [11]. Thus, 
for comparison with proton primary tm at energy 10^^ eV we must to shift the taken data to large 
dept at the least 71+0.2 ■ 85=87 (g/cm ^). To make visual comparison more convenient we shift 
the experimental data at 100 g/cm ^ in fig. 5. It is needed to have in view that both energy 
dependence and non-proton primaries wide the distribution significantly. So fig. 5 shows general 
agreement the simulated and experimental data. Remembering the descriptive ability of fen{x) 
(fig. 2) and quality data authors think that feni^) m^-Y be used to investigate distribution function 
of the depth of cascade maximum. 

The functions, P{tm), were obtained in above discussed way for all simulated cases primary 
protons: energies of 10^^ 10^° and lO^^ eV for zenith angles ^=0°, 24.6°, 44.4°, 60°. The full 
description of the distribution functions of tm is given by fig. 6 into compact form through mean 
values and components of variances. The mean values of tm are independent of angle, practically, 
and the elongation rate of proton fraction is seen clearly (^ 55g/cm'^), that is the same as [11]. The 
variance is decreasing at increasing the energy (^ 25% at all diapason). The position is natural 
and agrees with previous results for 45° [15]. Note, together with decreasing the variance one can 
see some increasing skewness of the distribution, but we think the conclusion may be unreliable 
(for detail see [5 (2007)]). Taking into account the quality of the data now, we can see scaling 
dependence of the distribution form at energy, only. The width of the band containing 50% of the 
probability is approximately 50g/cm^ (see also fig. 5). Thus, to define the primary energy, the 
parameter is inferior to s* = \/27rN^b\ 

Formal apparatus of distribution functions gives ability to obtain not only semi-quality estima- 
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tions (as above), but correct quantity estimations. To illustrate technical work with cumulative 
distribution presented by table 1 we will obtain the region of overlap the P{tm) for 6=A4A° at 
energy lO^^eV with one at energy lO^^eV. The parameters of the P{tm) take from left side of table 
2. The P{t^) at lO^o eV has left bound t„ > c - 2 ■ b'^/a = 807 - 2 ■ 8.5740 = 803A(g/cm'^). 
At the depth cumulative distribution (19) for 10^^ eV F(t^, 38.7, 22.1, 752.5) = 0.75. fl Thus, 
proton generated showers which energy differ more then in ten times can have the same depth of 
maximum with probability 0.25. To clear the problem, for example, one can create function of 
maximum probability at plane tm, - primary E, easily. 

6 Discussion 

As noted above, report [11] studies the cosmic ray composition in different energy ranges by com- 
paring the observed average tm with predictions from air shower simulations for different nuclei. 
The general situation is presented in figure 3 of the report. We added our results of distribution 
function of for proton generated showers: top straight line segment at 10^^ eV with arrows 
denoted ai, ar and lower conditional boundary for value of tm- See fig. 7. In our opinion, the 
difficulties in mass-spectra estimation due to not differences mean results of various models (as it 
was some years ago [15]), but wide band of physical fluctuations, i.e. wide width of distribution 
function. Therefore we hope the formal apparatus of distribution functions will be useful for the 
problem. 

At the highest energies, due to relict suppression of nuclear flux [1], it would be expected a 
significant proportion (10%-i- 50%) of the spectrum of cosmic rays would be photons [12]. The tm of 
photon generated shower is much greater then their nuclear counterparts. This is due to both causes: 
much lower multiplicity in particle production in electromagnetic - dominated photon showers than 
in the hadronic present for nuclear primaries and suppression of Bete-Heitler pair production by 
LPM effect [20], which is not important for other cosmic ray primaries. The influence of LPM 
effect on the form of cascade curve is well-known from electromagnetic cascade theory. Thus, both 
the results and tm can be used to discriminate between photon and nucleonic UHE primaries. 
More exactly estimation of the proportion may be carried out with help of formal apparatus of 
distribution functions (not only visual comparison of histograms). The increase of the statistic 
UHE events will give all needed for the ability. 

7 Summary 

The asymmetric form fen{x) = Aexp{—{x — cY /{a{x — c) + 26^)) inherits the convenience of inter- 
pretation of its parameters from exponential and normal distributions. The methods offer comple- 
mentary options to deal with incomplete and low-quality data, while remaining comfortable with 
the usual procedures for normal distribution. In particular, application of the presented apparatus 
for UHE-cascades enables to obtain the correct description of individual cascades generated by 
UHE protons and to investigate the distribution function of the depth of their maximum. 
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Table 1. The values of the functions for calculation of the parameters of fen{x, a, b) 
and for calculation of the values of Fcn{y- z). 



z = 


.lE-02 


.2E-02 


.5E-02 


.lE-01 


.2E-01 


.5E-01 


-lE+OO 


.2E+00 


.5E+00 


.lE+01 


A\a\ = 


.99900 


.99802 


.99509 


.99031 


.98109 


.95556 


.91826 


.85714 


.73233 


.61119 


x/a = 


.99950 


.99901 


.99757 


.99524 


.99081 


.97891 


.96232 


.93670 


.88952 


.84974 


a^a^ - 1 = 


.00000 


.00000 


.00001 


.00002 


.00010 


.00061 


.00235 


.00865 


.04304 


.12768 


of/'^f = 


.35914 


.35914 


.35914 


.35915 


.35918 


.35938 


.36006 


.36243 


.37406 


.39829 


-0.4 


.00000 


.00000 


.00000 


.00000 


.00000 


.00000 


.00000 


.00000 


.00000 


.00303 


-0.2 


.00000 


.00000 


.00000 


.00000 


.00000 


.00000 


.00000 


.00000 


.00602 


.07670 


0.0 


.00050 


.00099 


.00246 


.00485 


.00946 


.02221 


.04087 


.07143 


.13384 


.19441 


0.2 


.18168 


.18208 


.18326 


.18517 


.18880 


.19855 


.21229 


.23383 


.27550 


.31417 


0.4 


.33001 


.33034 


.33131 


.33287 


.33582 


.34367 


.35449 


.37084 


.39983 


.42329 


0.6 


.45146 


.45173 


.45252 


.45379 


.45620 


.46257 


.47125 


.48402 


.50486 


.51850 


0.8 


.55089 


.55111 


.55176 


.55280 


.55477 


.55996 


.56696 


.57707 


.59241 


.59981 


1.0 


.63230 


.63248 


.63301 


.63386 


.75565 


.63970 


.64538 


.65346 


.66493 


.66841 


1.2 


.69896 


.69910 


.69954 


.70023 


.70155 


.70500 


.70961 


.71610 


.72479 


.72586 


1.4 


.75353 


.75365 


.75400 


.75457 


.75599 


.75847 


.76222 


.76745 


.77409 


.77372 


1.6 


.79820 


.79830 


.79859 


.79906 


.79994 


.80224 


.80530 


.80953 


.81464 


.81346 


1.8 


.83478 


.83486 


.83510 


.83548 


.83620 


.83809 


.84058 


.84400 


.84797 


.84637 


2.0 


.86473 


.86480 


.86499 


.86531 


.86590 


.86744 


.86946 


.87224 


.87534 


.87357 


2.2 


.88925 


.88931 


.88947 


.88972 


.89020 


.89146 


.89312 


.89538 


.89780 


.89602 


2.6 


.92576 


.92580 


.92591 


.92608 


.92640 


.92724 


.92835 


.92984 


.93134 


.92978 


3.0 


.95024 


.95026 


.95033 


.95045 


.95067 


.95123 


.95197 


.95296 


.95390 


.95264 


3.5 


.96982 


.96983 


.96988 


.96995 


.97008 


.97042 


.97086 


.97146 


.97199 


.97110 


4.0 


.98169 


.98170 


.98173 


.98177 


.98185 


.98206 


.98232 


.98268 


.98299 


.98239 


5.0 


.99327 


.99327 


.99328 


.99329 


.99332 


.99340 


.99349 


.99362 


.99373 


.99347 


10.0 


.99995 


.99995 


.99996 


.99995 


.99996 


.99996 


.99995 


.99994 


.99996 


.99996 
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Continuation of tabic 1 



z 


.2E-I-01 


.5E-I-01 


.lE-f02 


.2E-f02 


.5E-I-02 


.lE-l-03 


.2E-t-03 


.5E-I-03 


.lE-l-04 


\ 11 = 


.48380 


.33318 


.24345 


.17518 


.11200 


.07949 


.05631 


.03566 


.02522 


x/a = 


.81431 


.78149 


.76708 


.75893 


.75368 


.75186 


.75093 


.75037 


.75021 


9/9 t 

CF^/a^ - 1 = 


.33690 


1.0448 


2.2749 


4.7632 


12.255 


24.752 


49.751 


124.75 


249.73 




.44472 


.54088 


.62683 


.70955 


.80079 


.85357 


.89369 


.93123 


.95083 


-50.0 


.00000 


.00000 


.00000 


.00000 


.00000 


.00000 


.00000 


.00000 


.00038 


-20.0 


.00000 


.00000 


.00000 


.00000 


.00000 


.00000 


.00061 


.02806 


.09296 


-10.0 


.00000 


.00000 


.00000 


.00000 


.00007 


.00968 


.05911 


.16908 


.25125 


-5.0 


.00000 


.00000 


.00000 


.00034 


.04221 


.12401 


.21285 


.30975 


.36337 


-4.0 


.00000 


.00000 


.00000 


.00551 


.08408 


.17493 


.25838 


.34257 


.38757 


-3.0 


.00000 


.00000 


.00043 


.03336 


.14743 


.23595 


.30795 


.37649 


.41218 


-2.5 


.00000 


.00000 


.00514 


.06431 


.18727 


.26980 


.33399 


.39378 


.42460 


-2.0 


.00000 


.00000 


.02503 


.10990 


.23207 


.30554 


.36072 


.41124 


.43708 


-1.6 


.00000 


.00218 


.06054 


.15706 


.27102 


.33526 


.38251 


.42532 


.44710 


-1.2 


.00000 


.02578 


.11696 


.21278 


.31222 


.36579 


.40459 


.43947 


.45715 


-1.0 


.00000 


.05387 


.15263 


.24335 


.33352 


.38130 


.41572 


.44656 


.46218 


-0.8 


.00060 


.09397 


.19256 


.27539 


.35517 


.39693 


.42689 


.45367 


.46722 


-0.6 


.02040 


.14462 


.23597 


.30861 


.37712 


.41267 


.43810 


.46079 


.47226 


-0.4 


.07842 


.20330 


.28199 


.34272 


.39929 


.42849 


.44933 


.46791 


.47730 


-0.2 


.16284 


.26715 


.32972 


.37742 


.42161 


.44436 


.46058 


.47504 


.48234 


0.0 


.25810 


.33341 


.37828 


.41241 


.44400 


.46025 


.47184 


.48217 


.48739 


0.2 


.35375 


.39971 


.42684 


.44740 


.46639 


.47615 


.48310 


.48930 


.49243 


0.4 


.44392 


.46418 


.47468 


.48212 


.48871 


.49202 


.49436 


.49643 


.49748 


0.6 


.52579 


.52545 


.52119 


.51632 


.51088 


.50784 


.50559 


.50356 


.50252 


0.8 


.59832 


.58261 


.56585 


.54978 


.53286 


.52358 


.51680 


.51067 


.50756 


1.0 


.66152 


.63511 


.60828 


.58230 


.55457 


.53923 


.52797 


.51778 


.51259 


1.2 


.71594 


.68272 


.64820 


.61371 


.57595 


.55475 


.53910 


.52488 


.51763 


1.4 


.76239 


.72544 


.68543 


.64388 


.59696 


.57013 


.55018 


.53196 


.52265 


1.6 


.80176 


.76342 


.71988 


.67268 


.61753 


.58535 


.56119 


.53903 


.52768 


1.8 


.83498 


.79692 


.75151 


.70005 


.63763 


.60038 


.57214 


.54608 


.53269 


2.0 


.86288 


.82627 


.78037 


.72591 


.65722 


.61521 


.58301 


.55311 


.53770 


2.2 


.88623 


.85184 


.80653 


.75025 


.67625 


.62982 


.59379 


.56011 


.54270 


2.6 


.92198 


.89309 


.85130 


.79429 


.71255 


.65830 


.61508 


.57405 


.55267 


3.0 


.94670 


.92356 


.88699 


.83228 


.74634 


.68571 


.63595 


.58788 


.56261 


3.5 


.96704 


.95026 


.92094 


.87179 


.78487 


.71830 


.66137 


.60498 


.57495 


4.0 


.97968 


.96794 


.94543 


.90333 


.81921 


.74889 


.68595 


.62184 


.58722 


5.0 


.99234 


.98697 


.97488 


.94706 


.87566 


.80370 


.73233 


.65476 


.61144 


10.0 


.99994 


.99988 


.99966 


.99848 


.98797 


.95804 


.89949 


.79746 


.72397 


20.0 


1.0000 


1.0000 


1.0000 


1.0000 


.99998 


.99949 


.99401 


.95324 


.88665 


50.0 


1.0000 


1.0000 


1.0000 


1.0000 


1.0000 


1.0000 


1.0000 


.99997 


.99858 
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Table 2. Results of method A and method (k=20) 
approximation of the probability density of the depth of cascade maximum for 9—AAA° 





A 


a 


b 


c 


P(A) 




A 


a 


b 


c 




10^9 


0.0147 


38.7 


22.1 


752.5 


0.88 


0.45 


0.0160 


41.5 


17.9 


749.5 


0.7 


1020 


0.0215 


40 


8.5 


807 


0.71 


0.25 


0.0202 


37.5 


12.2 


811 


0.23 


10^^ 


0.0255 


27 


11 


871 


0.72 


0.7 


0.0289 


27.8 


7.8 


869 


0.95 
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z=(2b/ar 



Fig. 1. Dependence of the ratio parts of variance (1 — crf/a^) on parameter of type of the distribu- 
tion. 
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Fig. 2. Example of standard forms of the distribution at a = 0, 5, 10; 6 = 1 and c = 0. The arrow 
indicates the mean value for each case. Here increasing a changes fen{x) from normal towards 
exponential. 
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Fig. 3. Interpretation of the parameters and approximation of a histogram by the hand method. 
For explanation of the construction see the text. 
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Fig. 4. Approximation of an individual shower by method B. Ordinates for N{t) (circles) are at 
left, ordinates for {t — tm)^ / \n{N{tm)/N{t)) (squares) are at right. 
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Fig. 5. Probability densities of the depth of UHE proton generated cascade maximum. Isolated 
fragment is representation of the A method result through mean value and components of variance. 
Experimental data is all-energy, all-angle and all-kind primary particle mixture, but is shifted on 
-1-100 g/cm^. For explanation see text. 
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Fig. 6. Distiibution functions of depth of simulated cascade maximum generated by protons. Logio 
energies (eV) are marked at the hues. Diameter points are accuracy of the mean values, arrows are 
parts of variance. Compare with fig. 5. 
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Fig. 7. Mean values of observed tm compared to predictions from hadronic interaction models. 
Experimental event numbers are indicated below each data point. Top strait line segment and 
arrows show distribution function of proton generated tm, lower segment is conditional bound of 
proton generated t^- 
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